Hydrodynamic interactions suppress deformation of suspension drops in Poiseuille 

flow 
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Evolution of a suspension drop entrained by Poiseuille flow is studied numerically at a low 
Reynolds number. A suspension drop is modelled by a cloud of many non-touching particles, ini- 
tially randomly distributed inside a spherical volume of a viscous fluid which is identical to the host 
fluid outside the drop. Evolution of particle positions and velocities is evaluated by the accurate 
multipole method corrected for lubrication, implemented in the hydromultipole numerical code. 
Deformation of the drop is shown to be smaller for a larger volume fraction. At high concentrations, 
hydrodynamic interactions between close particles signiflcantly decrease elongation of the suspension 
drop along the flow in comparison to the corresponding elongation of the pure-fluid drop. Owing 
to hydrodynamic interactions, the particles inside a dense-suspension drop tend to stay for a long 
time together in the central part of the drop; later on, small clusters occasionally separate out from 
the drop, and are stabilized by quasi-periodic orbits of the constituent non-touching particles. Both 
effects significantly reduces the drop spreading along the fiow. At large volume fractions, suspension 
drops destabilize by fragmentation, and at low volume fractions, by dispersing into single particles. 



I. INTRODUCTION 

Dynamics of a cloud of many non-touching micro- 
particles forming a suspension drop entrained by the 
Poiseuille flow in a microchannel is interesting for various 
industrial, biological and medical applications, such as 
e.g. transport in microfluidic devices [l| or narrow chan- 
nels , drug or gene delivery through the use of magnetic 
nano-particles Q or inhalator drug delivery 0, Q • The 
basic question, from a dynamical point of view, is how 
to control dispersion of the particles and the cloud defor- 
mation. This is important from the practical perspective, 
as in some aspects it is important that the particles or 
molecule clusters stay intact for a certain time or traveled 
distance, while in others fast dispersion is desired. 

It is known from other contexts 6] that hydrodynamic 
interactions which arise between non-touching particles 
moving in a viscous fluid under a low Reynolds number, 
in general tend to keep close particles together in groups. 
This phenomenon is well studied on theoretical, numeri- 
cal and experimental grounds 0,(1] for particle movement 
resulting from the gravitational field. In particular, an 
ensemble of initially randomly distributed particles on 
the average moves in the same way as a fluid drop set- 
tling down in a less dense host fluid. For a very long time, 
such a suspension drop behaves as a cohesive entity even 
though there is no surface tension to hold the suspended 
particles together. 

The clustering effect has been also observed for groups 
made of several particles only. In gravitational field, 
three particles can stay together for a very long time, 
what gives an indication to the evolution characteristics 
of larger particle clouds @ . In analogy, periodic motion 
of two particles under shear flow [To| forms a benchmark 
which we expect to influence dynamics of large clouds of 



particles entrained by an ambient fluid flow. In both con- 
texts, hydrodynamic interactions keep particles together 
for a long time. A similar mechanism has been observed 
in Ref. [11], where stability of particle clusters of an or- 
dered internal structure immersed in a shear flow was 
analyzed. 

A cloud of randomly distributed particles in an ambi- 
ent fluid flow, has yet not been investigated. This paper 
is therefore devoted to a study of the evolution of such 
a suspension drop entrained by a Poiseuille flow inside 
a two-wall channel, in the low-Reynolds-number regime. 
A suspension drop is modeled by a cloud of many non- 
touching particles, which are initially distributed ran- 
domly inside a given spherical volume of the fluid, iden- 
tical to the host fluid outside. In the course of evolution, 
the particles are free to move relative to each other and 
the surrounding channel walls. 

The goal is to analyze to what extent hydrodynamic in- 
teractions between the particles inside a suspension drop 
influence its evolution, and in particular, its elongation 
and dispersion along the flow. The strength of the hy- 
drodynamic interactions is tuned-up by the increasing 
volume fraction of the suspension. 

The structure of this paper is the following: section 
II contains the general system description followed by 
section III introducing briefly the numerical procedure 
used to calculate the hydrodynamic interactions. Section 
IV contains the main results and their discussion. The 
paper is summarized and concluded in section V. 



II. SYSTEM 

Consider a suspension drop made of N particles im- 
mersed in a viscous fluid identical to the fluid outside 
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the drop. Each individual particle is a hard sphere of di- 
ameter d. The stick boundary conditions on the particle 
surfaces are assumed. The particles cannot overlap, and 
do not interact with each other through electrostatic or 
magnetic forces or other direct interactions. The parti- 
cles are initially randomly distributed in a given spherical 
fluid volume of diameter D, with equal A^-particle prob- 
ability everywhere except overlapping. Volume fraction 
(f) of the particles is therefore given as 



The fluid is bounded by two infinite planar walls sepa- 
rated by a distance h, which is much larger than the drop 
radius. This configuration models a micro-channel of a 
width and length which are much larger than its height. 
Across the length of this channel a pressure gradient is 
exerted. In the absence of a drop, it would lead to the 
formation of a steady Poiseuille flow, 

vo = 4v™ z)//i^x. (2) 
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FIG. 1: Model system at initial moment of time. 

As illustrated in Fig. [1] the suspension drop is im- 
mersed in the ambient flow ([2]), which modifies the fluid 
velocity and pressure flelds, v(r) and p(r). A low- 
Reynolds-number flow is assumed and described by the 
steady Stokes equations [l^, : 

r/V^v - = 0, V • V = 0. (3) 

Therefore, translational velocities of the particle cen- 
ters, dvi/ dt, i = 1, N, are linear functions of the max- 
imal ambient flow velocity v^, 

dv • 

-^=v™C,(X), (4) 



where the coefficients Ci{X) are Cartesian vectors de- 
pending on the configuration of all the particle centers, 
X = (ri,r2, ...rjv). 



III. NUMERICAL PROCEDURE 

The Cartesian vectors Ci{X) are determined numeri- 
cally by the HYDROMULTIPOLE algorithm, which imple- 
ments the theoretical multipole method of calculating hy- 
drodynamic interactions between bodies [14i-,-17j within 
Stokesian dynamics. The defined cut-off parameter L for 
the multipole method was chosen L = 2 what means that 
24 multipole moments where calculated for each of the 
N particles. Such a choice of the multipole cut-off is suf- 
ficient to achieve a precision of the calculated velocities, 
normalized by v^, at least 5 x 10~^ when calculating ve- 
locities of particles in our system. This precision estimate 
was found by calculating the velocity of the suspension 
drop at i = for L ranging from 2 to 8 with cj) = 40% 
(dense drops, i.e. particles packed close together). 

The particle-wall interactions are incorporated by the 
single-wall superposition of hydrodynamic forces for the 
two- wall system |17| . This approximation is poor for nar- 
row channels when the initial diameter of the drop is 
comparable to the width of the channel but it performs 
well for wide channels Here we consider the case of 
relatively wide channels compared to the drop diameter 
when use of the single- wall superposition approximation 
is fully justified. 

The first analysis performed was the estimation of wall 
effects encountered. We compared the initial mean ve- 
locities of identical drops (same initial configuration of 
particles) in the same external Poiseuille fiow, using nu- 
merical HYDROMULTIPOLE procedures for an unbounded 
fluid and fluid between two parallel walls. At the initial 
moment of time, the difference between the two calcu- 
lated drop velocities, normalized by v^, was of the order 
of 2 X 10~^, i.e. smaller that the computational accuracy. 
Taking into account the above result, we used the nu- 
merical HYDROMULTIPOLE codcs for an unbounded fluid, 
and benefited from the three-time increase in calculation 
speed. 

Having calculated the instantaneous velocities of all 
particles, the evolution of the drop was determined by 
time stepping the set of coupled differential equations 
for each particle position. 

The number of independent simulation runs finally 
performed varied with the volume fraction. For (f> up to 
30%, ten independent initial configurations where con- 
sidered. For volume fractions 40% and 50% twenty sim- 
ulation runs where performed. This was motivated by 
larger fiuctuations observed for larger volume fractions. 
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IV. EVOLUTION OF A SUSPENSION DROP 

A. Model parameters 

The simulated model is described by three parameters, 
the channel height h, the number of particles in the sus- 
pension drop N and the volume fraction (p. The number 
of particles N in the drop was held constant, equal to 
A'' = 80. Our goal was to study how the drop evolution 
would change with the increased volume fraction, when 
the hydrodynamic interactions between the close parti- 
cles are enhanced. Increasing volume fractions, we also 
increased the size of the particles, d, according to Eq. ([T]) 
with a constant drop diameter D. At the same time, we 
kept a constant channel height h, so as the ratio D/h 
was also held constant. This ensured that the flow gradi- 
ent over a drop diameter was approximately equal for all 
volume fractions. Suspension drops of different volume 
fraction (i.e different diameter when composed of equal 
number of particles) could then be compared in terms of 
the influence the same flow had on the structure of the 
drop. We chose a low value D/h = 0.117, with the cor- 
responding values of h listed in Table HI This enabled us 
to make hydrodynamic interactions of the drop with the 
walls very weak, and therefore to focus on the hydrody- 
namic interactions between particles inside the drop. 



TABLE I: The channel height h/d used in our simulations for 
a given volume fraction (p, with TV = 80 and D/h = 0.117. 





h/d 


5% 


100.00 


10% 


79.37 


20% 


63.00 


30% 


55.04 


40% 


50.00 


50% 


46.42 



From now on, we normalize distances by the height h of 
the channel (i.e. the distance between the channel walls). 
Velocities are normalized by the maximum velocity of the 
flow, Vm, attained in the middle of the channel. Time t 
in the simulations is measured in h/vrn- All simulations 
were performed for times approximately up to i = 300. 



B. Initial moment of time 

The average velocity of a suspension drop 

M , N 



u 



M ^ N ■ 

j=i i=i 



(5) 



calculated at t = 0, is plotted in Fig. [2] as a function of 
volume fraction. Here A'^ = 80 is the number of particles 
in a drop and M = 100 is the number of independent sim- 
ulations performed for different random configurations of 



the particles, and [/^(j) is the velocity of particle i in sim- 
ulation j, divided by v„. Errors of the average veloc- 
ity(Eq. [5]) are smaller than the size of the plotted points. 
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FIG. 2: The dimensionless velocity of a suspension drop at 
t = averaged over 100 independent, random initial condi- 
tions (dots). Errors are smaller than the size of the dots. The 
dashed line is the limiting case of a volume of fluid having 
the same size and initial position as the considered suspen- 
sion drops. The dotted line corresponds to a suspension drop 
composed of non-interacting particles. 

In comparison, a spherical volume of pure fiuid identi- 
cal to the analyzed suspension drop moves at the initial 
moment of time with a velocity 



Uf 



l/jo(.).^^^(l-g), 



(6) 



where Q denotes the volume of a sphere of diameter D 
centered at (0,0,1/2) and vo(2;) — |vo|/v„i — 4z(l — z). 
This mean equals Uf = 0.997262 for D/h ^0.117 and is 
denoted by the dashed line in Fig. [2j 

The dotted curve is the velocity of a suspension drop 
composed of non-interacting particles, i.e. of a volume of 
fluid, corrected by the Faxen term, corresponding to the 
non-zero diameter d of particles. 



dU = 



24/i2 



ld2 
3/i2- 



(7) 



Therefore, according to Table HI the Faxen correction de- 
pends on volume fraction. For example, 5U = 1.3 x 10^'* 
for (j) = 0.40 (i.e. h = 50d). 

It can be noticed in Fig. [5] that a change of the suspen- 
sion volume fraction results in an overall small difference 
in mean velocity of the suspension drop, at most of the 
order of 10"'^. Nevertheless, the changes in evolution 
with the increase of volume fraction are substantial, as 
will be pointed out in the next section. 



C. Snap-shots from suspension drop evolution 

Although the difference in the mean velocity of the sus- 
pension drop changes only slightly with the change of 
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FIG. 3: Snap-shots from drop evolution simulations p^ - [23 |. All particles have been projected onto the y = plane. The red 
outline (given by Eq. HA3|) ) is the circumference of the projection of the instantaneous position of the fluid volume, which 
would be initially identical with the suspension drop, (enhanced online) 
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FIG. 4: Snap-shots from drop evolution simulations [UH^- AH particles have been projected onto the z = 1/2 plane. The red 
outline (given by Eqs. HA2IIA4|) ') is the circumference of the projection of the instantaneous position of the fluid volume, which 
would be initially identical with the suspension drop, (enhanced online) 
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its volume fraction, a quick look at the simulation snap- 
shots (Fig. El see [iliia] and Fig. H see dHH) shows 
a clear change in the behavior of all particles. As the 
volume fraction increases, particles tend to stay clustered 
together longer. The stretching of the initial shape of the 
suspension drop therefore highly depends of its volume 
fraction. In order to show this phenomenon clearly it is 
instructive to compare the positions of particles with the 
evolved shape of an initially spherical volume of fluid. 
The equation describing the circumferential surface of 
this fluid volume is 

(a;-Vo(z)i)2 = ^-z^-y\ (8) 
vo(z) = 4z(l-z). 

As simulations shown and discussed in this paper are 
presented in projections (?/ = or z = 1/2), the interior 
of the drop contained inside the above surface has to be 
accordingly projected. Therefore, 

1. The boundary of the fluid-drop projection on the 
y = plane is given by Eq. ([5]) with y = 0. 

2. The fluid-drop projection on the z = 1/2 plane is a 
superposition of circles shifted with respect to each 
other. Its boundary is therefore more sophisticated, 
see Eqs. (|Ml) - (|X4l) derived in Appendix E] 

The boundaries of the fluid-drop projections on the 
y — Q and z = 1/2 planes are plotted in Figs. [3] and 
m respectively, and compared with the snapshots of the 
underlying movies [click on figure to watch movie], pre- 
senting the suspension drops at the same time instant. 
The scale on the x, y and z axes is the same, to keep 
spherical shape of the particles and the initial volume of 
the drops. 

For small volume fractions the particles move in a simi- 
lar way as the pure fluid. Particles initially concentrated 
in a spherical drop get spread out evenly and the ini- 
tial shape of the drop is clearly stretched resembling the 
evolved shape of the corresponding spherical volume of 
pure fluid. No group formation is observed and hydrody- 
namic interaction between particles seems to be minor. 

As the volume fraction increases, the shape of the sus- 
pension drop at a given time instant is clearly less and 
less elongated by the flow than the evolved reference vol- 
ume of pure fluid, as shown in Figs. [3]|31 Notice that for 
the largest volume fraction (j) = 50%, the drop's shape is 
still quite close to a sphere. 

As the volume fraction increases, the evolution of the 
suspension drop changes significantly, and the particles 
have a tendency to stay clustered in a single group for a 
longer time, as illustrated in the movies linked to Figs.[3]- 
m For larger times, deformation of a dense suspension 
drop increases, and small groups of particles separate out 
from the main cluster. Particle groups which form in 
both tails, rotate relative to their center of mass, due to 
the gradient of the flow, in a similar way as two particles 
in a shear flow fl(5| . The average number of the particles 



in such a group increases with time and is larger for a 
higher volume fraction. 

D. Quantitative analysis 

By viewing the simulation results one can conclude 
that hydrodynamic interactions hold particles together 
more effectively, when these are packed in larger vol- 
ume fractions. This qualitative remark can be quantiflcd 
by analyzing the time-dependent dispersion of the parti- 
cle positions inside a drop, averaged over M simulations 
corresponding to different random initial configurations. 
The dispersion along x-direction is defined by the follow- 
ing formula, 




where A'^ = 80 is the number of particles in a drop 
and M is the number of simulations performed. The 
x-coordinates x^j-^ of each of the i — 1, . . . ^ N particles 
in each of the j = 1, . . . , M simulations are functions of 
time, and so is the dispersion itself. The same formula 
leads to the definition of ay and az once Xi(j) is exchanged 
for the appropriate coordinate. 

Fig. E] shows the dispersion ax of particle positions in 
the drop along the flow (x) direction, evolving in time, 
for different volume fractions. Color-coding, representing 
different volume fractions, has been chosen according to 
Table |lTl At a given time instant, the dispersion a^. is 
smaller for a larger volume fraction (j), and this effect is 
significant. 

0"x 




- — — — — — ^ — — — ^ — — — t 

50 100 150 200 



FIG. 5: Dispersion of particle positions in the drop along the 
flow (x) direction, evolving in time. Results averaged over all 
tested initial conditions. Color-coding see Table nil the lar ger 
volume fraction, the smaller is a^- The dashed line (</!> = 0%) 
corresponds to pure fluid, Eq. (|10p . and the dotted one to the 
non-interacting particles, Eq. (O- 

For comparison, the dashed line in Fig. El is the analyt- 
ical result calculated for a spherical volume of pure-fluid 
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TABLE II: Color-coding used to diflerentiate results for var- 
ious volume fractions 



color 





black 


5% 


red 


10% 


yellow 


20% 


green 


30% 


light blue 


40% 


dark blue 


50% 



evolving with the flow 
culated as 



0%). Its dispersion was cal- 



1 ^ im'^D^ 5L»2 

To' 



/l2 



where 



t I 



(10) 



(11) 



The dispersion ax of a pure-fluid drop is of course larger 
than the dispersion of a suspension drop. To investigate 
the effect of the excluded volume, the dispersion ax,ni 
of a suspension drop made of flctitious non-interacting 
particles is also evaluated, defined by Eq. ([9]) with the 
real particle positions x^jj{t) replaced by a;i(j),„i(i), the 
time-dependent positions of the fictitious particles, which 
would not interact hydrodynamically, and just translate 
along the Poiseuille fluid flow, vq^ with the Faxen velocity, 
Wo + SU, and the Faxen correction SU given by Eq. ([7]), 



Xt{j),m{t) = a;,(j)^„,(0) + {vo + SU)t. 



(12) 



The initial configurations of the fictitious particles are 
the same as those used in the suspension-drop evolution, 
Xi(j),ni(fi) = 3^i(j)(0)- The result is plotted in Fig. [5] with 
a dotted line. The curves corresponding to = 5%, 40% 
and 50% are practically superimposed. Comparison with 
the pure fluid and with the suspension drop indicates 
that elongation of the suspension drop is signiflcantly 
suppressed by strong hydrodynamic interactions between 
particles at larger volume fractions. 

In Fig. [S]a closeup of Fig. [5] for short times is shown. 
Further, the mean standard errors have been plotted 
(dotted lines) together with the original curves in order 
to show that the results obtained distinguish between dif- 
ferent volume fractions. For volume fractions 40% and 
50% the errors are smaller than the width of the curve. 

Systematically, for a larger volume fraction, the disper- 
sion ffx is smaller, if time is not too small, e.g. t > 5. No- 
tice that at t = the curves plotted in Fig. [HI do not over- 
lay each other. On first notice this might seem wrong. 
But when given further insight, the dispersion at t — 
is shown to change with the volume fraction, as plotted 
in Fig. [71 This graph shows the dispersion of particle 



0.08 



0.07 



0.06 



0.05 



0.04 



0.03 




FIG. 6: Dispersion of particle positions in the drop along the 
flow (x) direction evolving in time. Results averaged over all 
tested initial conditions. Top to botton: 5% to 50% respec- 
tively. Color-coding see Table Ull The dotted lines show the 
standard mean errors. For 40% and 50% the errors are smaller 
than the width of the curve. 

ff,(t=0) 



0.026 



0.0245 

0.024 

0.00001 



0=0 



0.0001 



0.001 
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FIG. 7: Dispersion of particle positions along the flow (x) 
direction calculated at f = 0. The dashed line is the limiting 
analytical solution. Error-bars are smaller then the size of the 
points. 



positions along the flow (x) direction calculated at t = 
for volume fraction ranging from 0.001% to 50%. The 
expected error of ax is equal the standard deviation of 
the dispersion within a single drop, divided by ^/M. For 
t = 0, the number of initial configurations randomly se- 
lected is 100000 and thus the error-bars are smaller then 
the size of the points. The non-monotonic dependence of 
the drop dispersion on the volume fraction is a strictly 
statistical effect due to excluded volume. In particular, 
notice that the dispersion for 5% is smaller than for 50% 
as also visible on Fig. [5] The dashed line is the limiting 
analytical solution at t — for a volume of pure fluid (i.e. 
= 0). Of course ay{t = 0) = cr^(t = 0) = axit = 0). 
We checked that for a larger number of particles in the 
suspension drop (and therefore for a smaller particle di- 
ameter d), the dispersion ax slightly increases. 

For larger times, the value of the drop dispersion in the 
flow direction decreases more than by a factor of 2, for 
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volume fractions changing form 5% to 50%. This gives 
a clear indication that hydrodynamic interactions tend 
to hold close particles together, and therefore suppress 
deformation of the drop along the flow. 

Dispersion of particle positions along the transverse 
(y) and (z) direction, respectively, evolving in time, is 
shown in the two last figures in this section - Fig. [S] and 
Fig. m Results have been averaged over all tested initial 
conditions. 



0.026 



0.024 



0.022 
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50% 
■40% 



FIG. 8: Dispersion of particle positions in the y direction, 
evolving in time. Results have been averaged over all tested 
initial conditions. The dashed line is the limiting result for 
a pure-fluid drop. The dotted lines show the standard mean 
errors. 



=0% 



0.026 




0.019 



FIG. 9: Dispersion of particle positions in the z direction, 
evolving in time. Results averaged over all tested initial con- 
ditions. The dashed line is the limiting result for a pure-fluid 
drop. The dotted lines show the standard mean errors. 



For comparison, we show also the corresponding dis- 
persion of the pure fluid. 



cr^^(i = 0) 0.02616. 



(13) 



In the transverse directions, the dispersion of the pure- 
fluid drop is constant in time. 

Note that in the direction transverse to the flow the 
drop gets contracted as it evolves. This effect is small. 



but clearly visible as the volume fraction of the drop in- 
creases. It is also caused by hydrodynamic interactions 
between close particles. 



E. Group formation 

When analyzing the simulation results at longer times, 
it becomes clear that for larger volume fractions, the 
particles which are lost from the main drop have a ten- 
dency to form small groups in the course of further evo- 
lution. This effect is readily visible in the movies linked 
to Figs. I3E1 and it is a clear indication of hydrodynamic 
interactions between close particles. 

In particular, for simulations with cj) > 30% groups of 
several particles where observed to stay together until the 
end of the simulation at t = 300. During this time each 
group rotated and particles interchanged places. This 
behavior of close particles resembles periodic trajectories 
of two particles in shear flow [10,]. 

We decided to check the group formation quantita- 
tively by introducing an algorithm working in the fol- 
lowing way: First, all particles are divided into two en- 
sembles, depending if their position z > 1/2 or 2; < 1/2. 
Then, within these two ensembles, positions of their cen- 
ters in the flow direction are compared by ordering. If 
between any two consecutive particle centers in such a 
sequence, the distance in the (x) flow direction is larger 
than a dimensional parameter g, then these particles are 
said to belong to two distinct groups. This grouping pa- 
rameter if fully arbitrary. We chose g = 3d. Different 
choices where checked - this discussion can be found in 
Appendix |B] 

Having in mind the definition of a group, the averaged 
time r of a first destabilization defined as formation of 
at least two groups of particles is evaluated and shown in 
Table IIIII and Fig. [TOl Clearly, for a larger volume frac- 
tion, hydrodynamic interactions between closer particles 
keep them together in a single cluster for a longer time. 



TABLE III: Averaged time r of flrst destabilization deflned 
as formation of at least two groups of particles. 



<^ 


T 


5% 


14.5 


10% 


20.5 


20% 


35.2 


30% 


44.0 


40% 


57.0 


50% 


90.0 



The error-bars in the figure correspond to the standard 
deviation of results calculated for the individual simula- 
tions. The line is a fit to the data 5% < < 40% given 
by the formula t = 9.0 + 120.3;/). The point for cj) = 50% 
is substantially above this fit, while all other results seem 
to be well in accordance with the linear behavior. 

The average number of distinct groups which are 
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0.1 0.2 0.3 0.4 0.5 

FIG. 10: Averaged time r of first destabilization as a func- 
tion of the drop volume fraction (j>. Error-bars are equal the 
standard deviation of the average value. The line is a linear 
fit to the data 5% < (/> < 40%. 



formed after time t = 50, 100, 150, 200 is listed in Ta- 
ble |IVl 



TABLE IV: Average number of distinct groups of particles 
after time t = 50, 100, 150, 200. 



<^ 


nt=50 


»T-t=100 


nt=i50 


?^t=200 


5% 


8.3 


18.1 


24.0 


30.0 


10% 


5.7 


13.3 


19.2 


22.0 


20% 


2.9 


9.0 


11.3 


13.3 


30% 


2.0 


5.8 


8.4 


10.1 


40% 


1.3 


4.2 


5.4 


6.7 


50% 


1.0 


2.1 


4.4 


6.2 



n 



5% 




, — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — ^ t 

50 100 150 200 

FIG. 11: Average number n of groups of particles as a function 
of time. Top-down: = 5%, 10%, 20%, 30%, 40% and 50%. 
The errors, calculated as the standard deviation of the mean, 
are smaller than the distance between the curves. 

Fig. [TT] shows the average number n of groups as a 
function of time, for all studied volume fractions. More 
dense drops split into a smaller number of groups. Taking 
into account the mean standard deviation as the error of 



these results, the dependence on the volume fraction is 
well established. Once again the difference between the 
case of 50% and 5% is substantial and suppresses 500%. 

The grouping phenomenon is therefore strictly corre- 
lated with the initial high concentration of particles. For 
large volume fractions, hydrodynamic interactions be- 
tween particles are strong and tend to cluster them. This 
is why more dense suspension drops destabilize slower. 
For a high volume fraction, particles stay in one group 
for a very long time, e.g. for times up to r = 90 if initially 
(j) = 50% (refer to Table |IV|. 

A higher volume fraction leads to a larger suspension 
viscosity. Therefore, the considered cloud of particles can 
be interpreted as a drop of a larger viscosity than a host 
fluid. In the absence of surface tension, the increase of 
the drop viscosity leads to its slower and smaller deforma- 
tion. This result gives an additional information to the 
numerical study of drop deformations at a finite capillary 
number, presented in Ref. (30| . 



V. CONCLUSIONS 

This paper was devoted to a numerical study of suspen- 
sion drop evolution in a Poiseuille flow of Stokesian fluid 
in a parallel-wall channel. The fluid inside the drop was 
the same as outside. The drop was initially centered on 
the axes of the channel, away from the walls. We studied 
the effect of the hydrodynamic interactions between sus- 
pended particles on the process of the drop deformation. 
Simulations where performed for a wide range of suspen- 
sion volume fractions, and compared to the evolution of 
a spherical volume of pure fluid, of the identical initial 
size and position in the Poiseuille flow as the suspension 
drop. 

The differences in evolution characteristics between the 
different volume fractions of suspension drops studied in 
this paper, are extensive. For a low volume fraction, the 
particles get dispersed evenly occupying an area covering 
approximately the evolved shape of a comparable vol- 
ume of pure fluid. Dense drops behave differently - for 
a long time they stay almost non-deformed; later, their 
evolution is dominated by formation of small groups of 
rotating particles, which are left behind the drop. The 
clustering is caused by hydrodynamic interactions. Rel- 
ative motion of close particles tends to hold them to- 
gether, pushing the flow out of the cluster and hindering 
its destabilization. The closer the particles are, the longer 
they stay together in a cluster. For example, at = 50%, 
the stretching of the initially compact drop, measured as 
the the dispersion of the particles along the flow, is two 
times smaller than at = 5%, and all the particles stay 
very close to each other in a compact group more than 6 
times longer. 

The results discussed in this paper have a clear impact 
on practical application. By changing volume fraction of 
suspension drops produced in microfluidic devices, or in- 
halators designed for drug delivery, one is able to control 



10 



the time-dependent dispersion of the particles entrained 
by the fluid flow. A large volume fraction will suppress 
the drop deformation and lead to clustering, while a small 
one to effective and even spreading of the particles. 
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Appendix B: Group detection 

The group detection algorithm, as described in sec- 
tion |IVE] depends on an arbitrary dimensional parameter 
g, which determines the inter-particle distance at which 
particles are said to belong to separate groups. We chose 
this parameter to be equal 3d. Correct conclusions can be 
drawn from the results of the algorithm, only if a change 
of g leads only to a quantitative not qualitative change 
of these results. The number of groups evaluated using 
the value oi g — 1.5d (i.e. half of the original value used 
in Fig. [TT|) is shown in Fig. [T^ 



Appendix A: Evolution of an initially spherical 
fluid-drop 

Volume of a pure-fluid drop remains the same dur- 
ing its evolution. The time-dependent circumferential 
surface of the fluid-drop has been specified in Eq. ([5]). 
Therefore, the boundary of the fluid-drop projection on 
the y — plane is given by Eq. ^ with y = 0, 



4/^2 



(Al) 



where yo{z) = 4z(l — z). An example of such a shape is 
shown in Fig. [31 

The circumference of the fluid-drop projection onto the 
plane z = 1/2 consists of the following two curves, 



XR.iy) = t + 



_D2_ 

4^ 



(A2) 



i.e. the right boundary of the drop projection, which is 
simply a half-circle, and 



XL{y) = niin |vo(z)i- 



(A3) 

i.e. the left boundary of the drop projection, which con- 
sists of two parts. 



XL{y) = ■ 



t- 



t- 



for y''>^-^ 



47? 



for y^<§^ 



1 

64F 



(A4) 



The curve trailing the area occupied by a projected, 
evolved shape of the spherical volume of fluid is there- 
fore in piece a circle given by Eq. (jA2p and the first 
part of Eq. (jA4[) . and piecewise a parabola given by the 
second part of Eq. (jA4l) . This curve and its derivative 
are continuous. An example of such a curve is shown in 

Fig.m 
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FIG. 12: Average number n of groups of particles as a function 
of time. Top-down: </> = 5%, 10%, 20%, 30%, 40% and 50%. 
The grouping parameter g — 1.5d. 

In Fig. [131 we have used g scaled with the width of 
the channel. In this case g has been chosen such that 
it is 3d for (p — 40% and for all other volume fractions 
it is scaled proportionally to the (dimensional) distance 
between the walls h, 




FIG. 13: Averaged number n of groups of particles as a func- 
tion of time, (j) — 5%, 10%, 50% as indicated, and top-down: 
= 20%, 30% and 40%. The grouping parameter g = 3ft/50. 



n-40% 50 



(Bl) 



11 



In both cases, we find that clearly drops of larger vol- eter g. Nevertheless, the general overall picture of the 
ume fraction tend to form fewer clusters. Of course, the evolution stays qualitatively the same, 
number of clusters depends on the choice of the param- 
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